#!/bin/bash
set -e

while getopts  ":db:p:" opts
do
        case  $opts  in
		b)
			asemmbly=$OPTARG
				;;
		d) anno_dbsnp=true;;
		p)              
			out_prefix=$OPTARG         
                ;;
		\?)
			echo `basename $0` [-d] in.vcf
			exit 65
			;;
        esac
done
shift $(($OPTIND - 1))

if [ -z $1 ]; then
	echo `basename $0` [-d] in.vcf
	exit 65
fi

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------

sample_name=$(grep -m1 '^#C' $1|awk '{print $NF}')

if test "$anno_dbsnp" = "true"; then
echo;echo; echo gatk VariantAnnotator dbsnp
java -Xmx$java_memory -jar $gatk \
-R $ref_genome \
-T VariantAnnotator \
-o $out_prefix.dbsnp.vcf \
-V $1 \
--dbsnp $data_path/ncbi/dbsnp/All_20150605.vcf.gz && mv $out_prefix.dbsnp.vcf $1
   
fi



format_af.sh $1
# format_af.sh -p2 1.indel.vcf 
anno_info_vcf.pl $1 $out_prefix.af.txt > $out_prefix.snp.af.vcf
# anno_info_vcf.pl 1.indel.vcf 2.af.txt > 1.indel.af.vcf

af_min=0.1
af_min_percent=`echo "$af_min * 100"|bc`
adp_min=6
adp_min_alt=2



echo 
echo 
echo gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.snp.filtered.vcf \
	--variant 1.snp.af.vcf \
	--filterExpression "FS > 60.0" \
	--filterName "FS_more_than_60" \
	--filterExpression "MQ < 40.0" \
	--filterName "MQ_less_than_40" \
	--filterExpression "HaplotypeScore > 200.0" \
	--filterName "HaplotypeScore_more_than_200" \
	--filterExpression "af < $af_min" \
	--filterName "Allele_frequency_less_than_$af_min_percent%" \
	--filterExpression "DP < $adp_min" \
	--filterName "Depth_less_than_$adp_min"
	
# --invertFilterExpression

:<<H
	
echo 
echo 
echo gatk VariantFiltration 
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.filtered.het.vcf \
	--variant 1.snp.filtered.vcf \
	--genotypeFilterExpression \
	--genotypeFilterName \
	--missingValuesInExpressionsShouldEvaluateAsFailing \


echo 
echo 
echo gatk SelectVariants 
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T SelectVariants \
	--variant 1.snp.filtered.vcf \
	-o 1.filtered.pass.vcf \
	-select 'vc.isNotFiltered()'
H


slt="vc.getGenotype(\"$sample_name\").getAD().1 > $adp_min_alt && vc.isNotFiltered()"


echo 
echo 
echo gatk SelectVariants
echo $slt
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T SelectVariants \
	--variant 1.snp.filtered.vcf \
	-o 1.filtered.pass.vcf \
	-select "$slt"


:<<!
	
# -select "FILTER == 'PASS'" # PASS means that the record is not filtered at all, so that expression won't work. You will need to use a more advanced JEXL expression (checking whether the VariantContext is filtered or not). See the docs of using JEXL expressions for more details.

het=\'vc.getGenotype\(\"$sample_name\"\).isHomRef\(\)\'
homvar=\'vc.getGenotype\(\"$sample_name\"\).isHomVar\(\)\'
echo 
echo 
echo gatk SelectVariants 
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T SelectVariants \
	--variant 1.filtered.pass.vcf \
	-o 1.pass.het.vcf \
	-select $het

echo 
echo 
echo gatk SelectVariants 
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T SelectVariants \
	--variant 1.filtered.pass.vcf \
	-o 1.pass.homvar.vcf \
	-select $homvar
	

# isHomRef


grep -v "#" my-favorite-snps.vcf | awk '{print $NF}' | awk -F ":" '{print $1}' | awk -F "/" '$1==$2 {print}'

cat 1.filtered.pass.vcf|grep -v "#" | awk '{if($1~/^[1-9]$/){print 0$1$2$4$5"/"$NF}else{print $1$2$4$5"/"$NF}}' | cut -d':' -f1|awk -F "/" '$2!=$3 {print}' >bulk.93.het.txt

!

# cat 1.filtered.pass.vcf|grep -v "#" | awk '{if($1~/^[1-9]$/){print 0$1$2"/"$NF}else{print $1$2"/"$NF}}' | cut -d':' -f1|awk -F "/" '$2!=$3 {print $1}' >sc.het.txt

# cat 1.filtered.pass.vcf|grep -v "#" | awk '{if($1~/^[1-9]$/){print 0$1$2"/"$NF}else{print $1$2"/"$NF}}' | cut -d':' -f1|awk -F "/" '$2==$3 {print $1}' >sc.hom.txt



cat 1.filtered.pass.vcf|grep -v "#" | awk '{if($1~/^[1-9]$/){print 0$1$2$3"/"$NF}else{print $1$2$3"/"$NF}}' | cut -d':' -f1|awk -F "/" '{if($2==$3){print $1 > "sc.hom.txt"}else{print $1 > "sc.het.txt"}}'
